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We show how Adjoint Algorithmic DifTerentiation (AAD) allows an extremely efficient calculation 
of correlation Risk of option prices computed with Monte Carlo simulations. A key point in the 
construction is the use of binning to simultaneously achieve computational efficiency and accurate 
confidence intervals. We illustrate the method for a copula-based Monte Carlo computation of 
claims written on a basket of underlying assets, and we test it numerically for Portfolio Default 
Options. For any number of underlying assets or names in a portfolio, the sensitivities of the option 
price with respect to all the pairwise correlations is obtained at a computational cost which is at 
most 4 times the cost of calculating the option value itself. For typical applications, this results in 
computational savings of several order of magnitudes with respect to standard methods. 
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One of the consequences of the current crisis of the Fi- 
nancial Markets is a renewed emphasis on rigorous Risk 
management practices. In order to quantify the financial 
exposure of financial firms, and to ensure an efficient cap- 
ital allocation, and more effective hedging practices, reg- 
ulators and senior management alike are insisting more 
and more on a thorough monitoring of Risk. Among all 
businesses, those dealing with complex, over the counter 
derivative securities are the ones receiving the most at- 
tention. 

A thorough calculation of the Risk exposure of port- 
folios of structured derivatives comes with a high oper- 
ational cost because of the large amount of computer 
power required. Indeed, highly time consuming Monte 
Carlo (MC) simulations are very often the only tool avail- 
able for pricing and hedging complex securities. Calcu- 
lating the Greeks, or price sensitivities, by 'Bumping' i.e., 
by perturbing in turn the underlying model parameters, 
repeating the simulation and forming finite difference ap- 
proximations results in a computational burden increas- 
ing linearly with the number of sensitivities computed. 
This easily becomes very significant when the models em- 
ployed depend on a large number of parameters, as it is 
typically the case. 

A particularly challenging task is the calculation of cor- 
relation Risk, i.e., the calculation of the sensitivites of a 
security with respect to some measure of the correlations 
among the random factors it depends on. Indeed, calcu- 
lating Risk with respect to all the independent pairwise 
correlations by Bumping requires repeating the MC sim- 
ulation a large number of times, e.g. increasing quadrat- 
ically with the number of random factors, and it is often 
unfeasible because of its high computational cost. 

Several alternative methods for the calculation of price 
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sensitivities have been proposed in the literature (for a 
review see e.g., [T]). Among these, the Pathwise Deriva- 
tive method [2] provides unbiased estimates at a compu- 
tational cost that may be smaller than the one of Bump- 
ing. However, in many problems the standard Pathwise 
Derivative method provides limited computational gains, 
especially when the contract priced has a complex pay- 
out [3 . A much more efficient implementation of the 
Pathwise Derivative method was proposed by Giles and 
Glasserman in Ref. [4 in the context of the Libor Mar- 
ket Model for European payouts, and recently general- 
ized to Bermudan options by Leclerc and co-workers [5]. 
These formulations express the calculation of the Path- 
wise Derivative estimator in terms of linear algebra op- 
erations, and utilize adjoint methods to reduce the com- 
putational complexity by rearranging appropriately the 
order of the calculations. 

Adjoint implementations can be seen as instances of 
a programming technique known as Adjoint Algorithmic 
Differentiation (AAD) 017]. In particular, as also dis- 
cussed in a forthcoming paper [8], AAD can be used as 
a design •paradigm to implement the Pathwise Derivative 
method, or the calculation of the sensitivities of any nu- 
merical algorithm, in full generality. In this paper we 
illustrate these ideas by discussing a specific application: 
the calculation of correlation Risk. We will begin by 
introducing the main ideas underlying Algorithmic Dif- 
ferentiation (AD), and the results on the computational 
efficiency of its two basic approaches: the Forward and 
Adjoint modes. 



Forward and Adjoint Algorithmic Differentiation 

Both the Forward and Adjoint mode of AD aim at 
calculating the derivatives of a computer implemented 
function. They differ by the direction of propagation of 
the chain rule through the composition of instructions 
representing the function. To illustrate this point, sup- 
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pose we begin with a single input a, and produce a single 
output z after proceeding through a sequence of steps: 



The Forward (or Tangent) mode of AD (FAD) defines 
ii to be the sensitivity of u to changes in a, i.e., 

. _ du 
da 

If the intermediate variables u and v are vectors, v is 
calculated by differentiating the dependence of u on u so 
that 



Vi 



dvi 



^ du., 



Applying this to each step in the calculation, working 
from left to right, we end up computing i, the sensitivity 
of the output to changes in the input. Note that if we 
have more than one input, then we need to calculate the 
sensitivity to each one in turn, and so the cost is linear 
in the number of input variables. 

Instead, the Adjoint (or Backward) mode of AD 
(AAD) works from right to left. Using the standard AD 
notation, u is defined to be the sensitivity of the output 
z to changes in the intermediate variable u, i.e. 



dz 
dui 



Using the chain rule we get. 



dz 
du. 
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which corresponds to the adjoint mode equation 



E 



dvj_ 
du,. 



Starting from z = 1, we can apply this to each step in 
the calculation, working from right to left, until we ob- 
tain a, the sensitivity of the output to each of the input 
variables. 

In the Adjoint mode, the cost does not increase with 
the number of inputs, but if there is more than one output 
then the sensitivities for each output have to considered 
one at a time and so the cost is linear in the number 
of outputs. Furthermore, because the partial derivatives 
depend on the values of the intermediate variables, one 
first has to compute the original calculation storing the 
values of all of the intermediate variables such as u and 
w, before performing the Adjoint mode sensitivity calcu- 
lation. 

In the above description, each step can be a distinct 
high-level function, or specific mathematical operations, 
or even an individual instruction in a computer code. 



This last viewpoint is the one taken by computer scien- 
tists who have developed tools which take as an input a 
computer code to perform some high-level function, 

V = FUNCTION([/) 

and produce new routines which will either perform the 
standard sensitivity analysis 

V = FUNCTION_D(C/, U) 

with suffix _D for "dot" , or its adjoint counterpart 

U = FUNCTI0N_B(;7, V) 

with suffix _B for "bar" [Tl] . 

One particularly important theoretical result is that 
the number of arithmetic operations in the adjoint rou- 
tine FUNCTION_B is at most a factor 4 greater than in 
FUNCTION [6]. As a result, it is possible to show that 
the execution time of FUNCTION_B is bounded by approx- 
imatively 4 times the cost of execution of the original 
function FUNCTION. Thus, one can obtain the sensitivity 
of a single output to an unlimited number of inputs for 
little more work than the original computation. 

While the application of such automatic AD tools to 
large inhomogeneous pricing softwares is challenging, the 
principles of AD can be used as a programming paradigm 
that can be used to design the Forward or Adjoint of 
any algorithm (possibly using automatic AD tools for the 
implementation of smaller, simpler components). This is 
especially useful for the most common situations where 
pricing codes use a variety of libraries written in different 
languages, possibly linked dynamically. These ideas will 
be discussed at length in Ref. [S]. 



AAD and the Pathwise Derivative method for 
Correlation Risk 

In this paper, we consider options pricing problems 
that can be expressed as an expectation value of the form 



y = Eq P{X) , 



(1) 



where X — (Xi, . . . , X^Y represents the state vector of 
N market factors (e.g., stock prices, interest rates, foreign 
exchange pairs, default times etc.), P{X) is the (possibly 
discounted) payout function of a security contingent on 
their future realization, and Q — Q(A) represents a risk 
neutral probability distribution [9 according to which 
the components of X are distributed. Although the pro- 
posed method easily generalizes to other kinds of joint 
distributions, here we consider a A'^-dimensional Gaus- 
sian copula as a model for the co-dependence between 
the components of the state vector, namely a joint cu- 
mulative density function of the form 



= $Ar($-l(Mi(Ai)), . . . , $-l(Afw(XAr));p) (2) 
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where ^n{Zi, . . . , Zn; p) is a A^-dimensional muhivari- 
ate Gaussian distribution with zero mean, and a N x N 
positive semidefinite correlation matrix p; is the 

inverse of the standard normal cumulative distribution, 
and Mi{Xi), i = 1, . . . , N, are the Marginal distributions 
of the underlying factors, typically implied from the mar- 
ket prices of liquid securities. 

The expectation value in ([TJ can be estimated by 
means of MC by sampling a number Nmc of random 
replicas of the underlying state vector . . . , X[AfMc]j 
according to the distribution Q(X), and evaluating the 
payout P{X) for each of them. This leads to the central 
limit theorem [ID] estimate of the option value V as 



The Pathwise Derivative method allows the calculation 
of the sensitivities of the option price V Q with respect 
to a set of Ng parameter 9 — {9i, . . . , OnbTi say 



dV{0) 
d6k 



_d_ 



P{X) 
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by defining appropriate estimators, say 6}^{X[iMc\), that 
can be sampled simultaneously in a single MC simula- 
tion. This can be achieved by observing that whenever 
the payout function is regular enough (e.g., Lipschitz- 
continuous, see Ref. [T]), and the distribution P(Z') does 
not depend on 0, one can rewrite Eq. ([s]) by taking the 
derivative inside the expectation value, as 



A'mc 



V ■ 



MC 



(3) 



■iMC = l 



dV{9) 
d9k 



dP {X) 
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with standard error S/ %/ A^mc , where — Eiq[P {X) ] — 
Eiq[P {X)]'^ is the variance of the sampled payout. 

In the Gaussian model above the dependence be- 
tween the underlying factors is determined by the 
correlation of a set of jointly normal random vari- 
ables Z — {Zi,...,ZffY distributed according to 
4>Ar(Zi, . . . , Zn; p). Each Zi is distributed according to a 
standard normal distribution so that '^{Zi) is a uniform 
random variable in (0,1) and Xi — M~^($(Zi)) is dis- 
tributed according to Mi. The sampling of the A'' jointly 
normal random variables (Zi, . . . , Z^r) is efficiently im- 
plemented by means of a Cholesky factorization of the 
correlation matrix. The Cholesky factorization produces 
a lower triangular NxN matrix C such that p — CC^ so 
that one can write Z — CZ where Z — (Zi, . . . , Zj^Y is 



The calculation of Eq. ^ can be performed by apply- 
ing the chain rule, and computing the average value of 
the so-called Pathwise Derivative estimator 



dPjX) 
d9k 



N 



sr^ dP{X) ^ dX, 
^ dX, d9k 
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The standard pathwise implementation corresponds to 
a Forward mode sensitivity analysis. Applied to steps 1-4 
(since the normal variates Z do not depend on any input 
parameters), this gives for each sensitivity: 

Step If Calculate Z = C Z where C is the sensitivity of 
C with respect to a given entry of the correlation 
matrix. 



(Zi) Zi, 



1, 



,Ar. 



a A^ dimensional vector of independent standard normal Step 2f Set f/j 
random variables. These observations naturally translate 

into the standard algorithm to generate MC samples of Step 3f Set Xi — Ui / mi(Xi), i — 1, . . . ,N. 
X according to ([2|, namely 



N 



Step Generate a sample of A^ independent standard nor- 
mal variates, Z = (Zi, . . . , Z^Y . 

Step 1 Correlate the components of Z by performing the 
matrix vector product Z — CZ. 

Step 2 Set Ui = $(Z,), i = 1, . . . , A. 

Step 3 Set Xi = M^^{Ui), i = l,...,N. 



BP 

Step 4f Calculate P = V — X 



Step 4 Compute the payout estimator P{Xi 



,Xn)- 



Correlation Risk can be obtained in an highly efficient 
way by implementing the so-called Pathwise Derivative 
method according to the principles of AAD [3], IH* . It is 
convenient to first express the expectation value as being 
over P(Z), the distribution of independent Z used in the 
MC simulation, so that 



(4) 



y = Eq p{x) =Ep Pixiz) 



The point of this subtle change is that P{Z) does not 
depend on the correlation matrix p, whereas Q(A) does. 



, , dXi 

Here (j){x) = d^{x)/dx is the standard normal proba- 
bility density function, and mi{x) = dMi{x) / dx is the 
probability density function associated with the marginal 
Mi{x) of the i-ih random factor. 

As anticipated, the computational cost of the Forward 
Pathwise Derivative method scales linearly with the num- 
ber of sensitivities computed Ng, i.e., the same scal- 
ing of finite difference approximations of the derivatives 
dg^EQ[P{X)]. As a result in many situations, typically 
involving complex payouts, the standard implementation 
of the Pathwise Derivative method offers a limited com- 
putational advantage with respect to Bumping 

In contrast, AAD allows in general a much more ef- 
ficient implementation of the Pathwise Derivative esti- 
mators ([7). Indeed, as an immediate consequence of the 
computational complexity results introduced in the pre- 
vious Section, it can be shown [8] that AAD allows the 
simultaneous calculation of the Pathwise Derivative esti- 
mators for any number of sensitivities at a computational 



4 



Chalesky 


ta(/*In*/ rho,/*In*/ L b,/*Out*/ rho b) 




1 

// 


n = Dimension of Correlation Matrix rho 




// 


Forward Sweep 




for 


(unsigned long i = 0; i < n; i++) { 






for (unsigned long j = i; j < n; j++) [ 






sumd, j) - rho(i, j) ; 






for (long k = i - 1 ; k >= D ; k — ) 






sum(i, j) — L(i, k) • l(j. 


k) ; 




if (i " jl ( 






1(1, i) - sgrt(sum(i. 


1) ) ; 




) else ( 






L(j, i) - sumd, j) / 

) 

) 


1(1, il; 


) 

II 


> 

Backward Sweep 




for 


(long i . n - 1; 1 >. 0; i— ) ( 






for (long j = n-1; j >= 1; j--) { 






double sum b = 0.0; 






if (i — j) 1 






if (sum(l, j) — 0.0) { 






sum b = 0.0; 






) else ( 






sum b=Lb(i, j)/( 2 


* L ( 1 , j ) ) ; 




) 

L b(i, j) = 0.0; 


1 




1 else \ 


1 




sum b = L b(j, i)/L(i, i) ; 






L b(i, 1) — surnd, jl • su 


i_b / L(i, 1) ; 




L b(j, i) =0.0; 






1 

for (long k = i - 1; k >= 0; k — ) 






L b(i, k) -= L(j, k}*sum b; 






L b(j ,k) -= L(i, k)*sum b; 






1 

rho b(i, j) += sum b; } 

1 




) 







given by 



FIG. 1: Adjoint of the Cholesky factorization. The Forward 
sweep is an exact replica of the original factorization. 



cost which is a small multiple (of order 4) of the cost of 
evaluating the original payout estimator. As a result, one 
can calculate the MC expectation of an arbitrarily large 
number of sensitivities at a small fixed cost. 

Although AAD can be applied for virtually any model 
and payout function of interest in Computational Fi- 
nance - including path-dependent and Bermudan options 
- here we will concentrate on the calculation of correla- 
tion sensitivities in a Gaussian copula framework. In 
general, for the reasons mentioned in the previous Sec- 
tion, the AAD implementation of the Pathwise derivative 
method contains a forward sweep - reproducing the steps 
followed in the calculation of the estimator of the option 
value P{X) - and a backward sweep. As a result, the ad- 
joint algorithm consists of adjoint counterparts for each 
of the Steps 1-4 above executed in reverse order, plus the 
adjoint of the Cholesky factorization. 

The first step consists in the evaluation of the adjoint 
of step 4 of the Forward sweep, calculating the derivatives 
of the Payout with respect to the components of the state 
vector 



dPjX) 
dX, 



(8) 



with k — 1, . . . ,N. These derivatives can be calculated 
efficiently using AAD, as discussed in Ref. [3]. 

In turn, the adjoint of Step 3 of the Forward sweep is 



Uk^Mk\Uk,Xk) 



Xh 



mk{Xk) 



(9) 



for k — 1,...,N. The vector U is then mapped into 
the adjoint of the correlated standard normal variables 
Z through the counterpart of Step 2 



Zk = HZk,Uk)^Uk^{Zk) 



(10) 



The adjoint of Step 1 performing the matrix vector prod- 
uct Z = CZ reads 



N 



a. 



-^^^ — Z}^ — Zj 



fc=l 



(11) 



or C = ZZ^ . By applying the chain rule, it is straight- 
forward to realize that the adjoint w of each intermedi- 
ate variable w in the succession of Steps 0-4 represents 
the derivative of the Payout estimator with respect to w, 
or w = dP/dw. In particular the quantities Cij calcu- 
lated at the end of the adjoint of Step 1 represent the 
derivatives of the payout estimator with respect to the 
the entries of the triangular Cholesky matrix C, namely 
the pathwise estimator ([T]) with 6^ = Cij. 

In summary, the AAD implementation of the Pathwise 
Derivative Estimator consists of Step 1-4 described above 
(forward sweep) plus the following steps of the backward 
sweep: 

Step 5 Evaluate the Payout adjoint X/. = dP/dX^, for 
fc = 1 N. 



Step 6 Calculate Uk = Xk/mk{M^^{Uk)), fc = 1, . 



.,N. 



Step 7 Calculate Zk = Uk(f>{Zk), k = I, 
Step 8 Calculate C = ZZK 



At this point in the calculation, there is an interesting 
complication. The natural AAD approach would average 
the values of C from each of the MC paths. This average 
C can be converted into derivatives with respect to the 
entries of the correlation matrix p by means of the adjoint 
of the Cholesky factorization jTT], namely a function of 
the form 



p = CHOLESKY _B(p, C) 



providing 



N 



/,m— 1 '-^ 



(12) 



(13) 



The pseudocode for the adjoint Cholesky factorization 
is given in Fig. [l] By inspecting the structure of the 
pseudocode it appears clear that its computational cost is 
just a small multiple (of order 2) of the cost of evaluating 
the original factorization. Indeed, the adjoint algorithm 
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essentially contains the original Cholesky factorization 
plus a backward sweep with the same complexity and a 
similar number of operations. 

The complication with this implementation is that it 
gives an estimate for the correlation risk, but it does not 
provide a corresponding confidence interval. An alterna- 
tive approach would be to convert C to p for each indi- 
vidual path, and then compute the average and standard 
deviation of p in the usual way. However, the numerical 
results will show that this is rather costly. An excellent 
compromise between these two extremes is to divide the 
Nmc paths into Nf, 'bins' of equal size. For each bin, 
an average value of C is computed and converted into a 
corresponding value for p. These Nb estimates for p can 
then be combined in the usual way to form an overall 
estimate and confidence interval for the correlation risk. 

The computational benefits can be understood by con- 
sidering the computational costs for both the standard 
evaluation and the adjoint Pathwise Derivative calcula- 
tion. In the standard evaluation, the cost of the Cholesky 
factorization is 0{N^), and the cost of the MC sampling 
is 0{NmcN'^), so the total cost is 0{N^ + NmcN^)- 
Since N^c is always much greater than N, the cost of the 
Cholesky factorization is usually negligible. The cost of 
the adjoint steps in the MC sampling is also 0{NmcN^), 
and when using Nt bins the cost of the adjoint Cholesky 
factorization is 0{Ni,N^). To obtain an accurate con- 
fidence interval, but with the cost of the Cholesky fac- 
torisation being negligible, requires that Ni, is chosen so 
that 1 < A^b < Nmc/N. Without binning, i.e., us- 
ing Nf, — Nmci the cost to calculate the average of the 
estimators (13) is 0{Ni\icN^), and so the relative cost 
compared to the evaluation of the option value is 0{N). 

The binning procedure described above can be gener- 
alized to any situation in which the standard solution 
procedure involves a common preprocessing step before 
any of the path calculations are performed. Other ex- 
amples would include calibration of model parameters to 
market prices, or a cubic spline construction of a local 
volatility surface. In each case, there is a linear relation- 
ship between the forward mode sensitivities before and 
after the preprocessing step, and therefore a linear rela- 
tionship between the corresponding adjoint sensitivities. 

The algorithm described above can be applied when- 
ever the option pricing problem can be formulated as an 
expectation value over a set of random factors whose dis- 
tribution is modelled as a Gaussian copula. This include 
in general a variety of Basket Options common across all 
asset classes, or structured swaps whose coupon depends 
on a specific observation of a set of correlated rates. In 
addition, the same ideas can be extended to the simula- 
tion of correlated diffusion processes [S]. 



Numerical Tests 

As a numerical test ground we consider the case of 
Basket Default Options 12J. In this context, the random 
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FIG. 2: Ratios of the CPU time required for the calculation 
of the option value, and correlation Greeks, and the CPU time 
spent for the computation of the value alone, as functions of 
the number of names in the basket, for Nmc = 10^. Sym- 
bols: Bumping (one-sided finite differences) (triangles), AAD 
without binning (i.e. A^f, = Nmc) (stars), AAD with binning 
{Ni, = 20) (empty circles). Lines are guides for the eye, and 
the MC uncertainties are smaller than the symbol sizes. 



factors Xi represent the default time of the i-th name, 
e.g., the time a specific company in a reference pool of N 
names fails to pay one of its liabilities as specified by the 
terms of the contract priced. In particular, in a n-th to 
default Basket Default Swap one party (protection buyer) 
makes regular payments to a counterparty (protection 
seller) at time Ti , . . . , Tm < T provided that less than 
n defaults events among the components of the basket 
are observed before time Tm- On the other hand, if n 
defaults occur before time T, the regular payments cease 
and the protection seller makes a payment to the buyer 
of (1 — i?i) per unit notional, where Ri is the normalized 
recovery rate of the i-th asset. The value at time zero 
of the Basket Default Swap on a given realization of the 
default times ri, . . . , rjv, i.e., the Payout function, can be 
therefore expressed as 

P(ti, . . . , Tat) = Pprotiri, . . . ,Tn) — Pprern{Ti, . . . , Tn) 

(14) 

i.e., as the difference between the so-called protection and 
premium legs. The value leg is given by 

PprotiTi, . . . ,TAr) = (1 - i?„)I?(r)I(r < T) , (15) 

where i?„ and r are the recovery rate and default time 
of the n-th to default, respectively, D{t) is the discount 
factor for the interval [0, t] (here we assume for simplicity 
uncorrelated default times and interest rates), and I(t < 
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T) is the indicator function of the event that the rt-th 
default occurs before T. The premium leg reads instead, 
neglecting for simplicity any accrued payment, 



p 



,(ri, . . . ,tn) 



L(t) 

E 



SkDiTk) 



(16) 



where L{t) = max[A; G {1, . . . , Af}/Tfc < r], and Sk is 
the premium payment (per unit notional) at time T^.. 
In order to apply the Pathwise Derivative method to 



the payout above, the indicator functions in ( 16 1 and 
(15), need to be regularized [U [12]. One simple and 



practical way of doing that is to replace the indicator 
functions with their smoothed counterpart, at the price 
of introducing a small amount of bias in the Greek esti- 
mators. For the problem at hand, as it is also generally 
the case, such bias can be easily reduced to be smaller 
than the statistical errors that can be obtained for any 
realistic number of MC iteration Nmc (for a more com- 
plete discussion of the topic of payout regularization see 
Refs. [3, 8, 13 ). 

The remarkable computational efficiency of AAD is il- 
lustrated in Fig. |2] for the Second to Default Swap. Here 
we plot the ratio of the CPU time required for the cal- 
culation of the value of the option, and all its pairwise 
correlation sensitivities, and the CPU time spent for the 
computation of the value alone, as functions of the num- 
ber of names in the basket. As expected, for standard 
finite-difference estimators, such ratio increases quadrat- 
ically with the number of names in the basket. Already 
for medium sized basket {N ~ 20) the cost associated 
with Bumping is over 100 times more expensive than the 
one of AAD. 

Nevertheless, at a closer look (see the inset of Fig. [2]), 
the relative cost of AAD without binning is 0{N), for 
the reasons explained earlier. However, when using Ni, = 
20 bins the cost of the adjoint Cholesky computation is 



negligible and the numerical results show that all the 
Correlation Greeks can be obtained with a mere 70% 
overhead compared to the calculation of the value of the 
option. This results in over 2 orders of magnitude savings 
in computational time for a basket of over 40 Names. 

Conclusions 

In conclusion, we have shown how Adjoint Algorith- 
mic Differentiation allows an extremely efficient calcula- 
tion of correlation Risk in Monte Carlo. The proposed 
method relies on using the Adjoint mode of Algorithmic 
Differentiation to organize the calculation of the Pathwise 
Derivative estimator, and to implement the adjoint coun- 
terpart of the Cholesky factorization. For any number of 
underlying assets or names in a portfolio, the proposed 
method allows the calculation of the complete pairwise 
correlation Risk at a computational cost which is at most 
4 times the cost of calculating the option value itself, re- 
sulting in remarkable computational savings with respect 
to Bumping. We illustrated the method for a Gaussian 
copula-based Monte Carlo computation, and we tested 
it numerically for Portfolio Default Options. In this ap- 
plication, the proposed method is 100 times faster than 
Bumping for 20 names, and over 1000 times for 40 names. 
The method generalizes immediately to other kind of El- 
liptic copulas, and to a general diffusive setting. In fact, 
it is a specific instance of a general AAD approach to the 
implementation of the Pathwise Derivative method that 
will be discussed in a forthcoming publication |8j. 
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